---
title: "BeckHall 2018 Meta-Analysis Code"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(metafor)
library(forestplot)
library(ggpubr)
library(ggplot2)

```

##Histograms from Overall Meta-analysis

```{r}

MetaN=read.csv("BeckHall_2018_OverallMetaAnalysis_NData.csv")
MetaP=read.csv("BeckHall_2018_OverallMetaAnalysis_PData.csv")
MetaNP=read.csv("BeckHall_2018_OverallMetaAnalysis_NPData.csv")

ggsave("Histograms_plot.png")

Nplot=ggplot(data=MetaN, aes(MetaN$N_LRR)) + 
  geom_histogram(breaks=seq(-7, 6, by = .25), 
                 col="black", 
                 fill="grey", 
                 alpha = .2) + 
  labs(x="N Log Response Ratio", y="Frequency")+
  theme(panel.grid.major = element_blank(), panel.grid.minor=element_blank(), axis.title.x=element_text(size=12), axis.title.y=element_text(size=12), axis.text.y=element_text(colour="black"), axis.text.x=element_text(colour="black"))

Pplot=ggplot(data=MetaP, aes(MetaP$P_LRR)) + 
  geom_histogram(breaks=seq(-3.5, 4, by = .25), 
                 col="black", 
                 fill="grey", 
                 alpha = .2) + 
  labs(x="P Log Response Ratio", y="Frequency")+
  theme(panel.grid.major = element_blank(), panel.grid.minor=element_blank(), axis.title.x=element_text(size=12), axis.title.y=element_text(size=12), axis.text.y=element_text(colour="black"), axis.text.x=element_text(colour="black"))


NPplot=ggplot(data=MetaNP, aes(MetaNP$NP_LRR)) + 
  geom_histogram(breaks=seq(-6, 6, by = .25), 
                 col="black", 
                 fill="grey", 
                 alpha = .2) + 
  labs(x="NP Log Response Ratio", y="Frequency")+
  theme(panel.grid.major = element_blank(), panel.grid.minor=element_blank(), axis.title.x=element_text(size=12), axis.title.y=element_text(size=12), axis.text.y=element_text(colour="black"), axis.text.x=element_text(colour="black"))

figure <- ggarrange(Nplot, Pplot, NPplot, labels = c("A", "B","C"),
                    ncol = 2, nrow = 2)
annotate_figure(figure,
                top = text_grob("Histograms of Effect Sizes", face = "bold", size = 18))

```

##Microbial Meta-Analysis Chlorophyll a Effect Sizes

```{r }

Data=read.csv("BeckHall_2018_MicrobeMetaAnalysis_LogResponseRatios.csv")
NEffect=subset(Data, Treatment=="N")
PEffect=subset(Data, Treatment=="P")
NPEffect=subset(Data, Treatment=="NP")


##Multi-level Regressions with Experiments Nested within Stream for Chla (Overall Effect Sizes, Data 1)

Chla_N <- rma.mv(Chla_LRR, Chla_LRR_Var, random = ~ 1 | Site, data=NEffect)
summary(Chla_N)

Chla_P <- rma.mv(Chla_LRR, Chla_LRR_Var, random = ~ 1 | Site, data=PEffect)
summary(Chla_P)

Chla_NP <- rma.mv(Chla_LRR, Chla_LRR_Var, random = ~ 1 | Site, data=NPEffect)
summary(Chla_NP)

#Forestplot

tabletext <-cbind (
  c("Nitrogen (n=51)", "Phosphorus (n=45)", "Nitrogen + Phosphorus (n=41)"))

forestplot(labeltext=tabletext, 
           mean=c(0.5018, 0.4426, 1.0391), 
           lower=c(0.2357, 0.1173, 0.6489), upper=c(0.7678, 0.7678, 1.4293),
           title="Chla Effect Sizes by Nutrient Treatment",
           xlab="LRR",
           xticks=c(0, 0.25, 0.5, 0.75, 1, 1.25, 1.5),
           txt_gp=fpTxtGp(label=gpar(cex=1.25),
                          ticks=gpar(cex=1.25),
                          xlab=gpar(cex=1.25),
                          title=gpar(cex=1.5)),
           col=fpColors(box="black", lines="black", zero = "gray50"),
           cex=0.5, lineheight = "auto", boxsize=0.15, colgap=unit(6,"mm"),
           lwd.ci=1, ci.vertices=TRUE, ci.vertices.height = 0.05)


```

##Microbial Meta-Analysis AFDM Effect Sizes

```{r}
##Multi-level Regressions with Experiments Nested within Stream for AFDM (Overall Effect Sizes)

NEffect2=subset(NEffect, AFDM_LRR_Var<2) #removing outliers
AFDM_N <- rma.mv(AFDM_LRR, AFDM_LRR_Var, random = ~ 1 | Site, data=NEffect2)
summary(AFDM_N)

AFDM_P <- rma.mv(AFDM_LRR, AFDM_LRR_Var, random = ~ 1 | Site, data=PEffect)
summary(AFDM_P)

NPEffect2=subset(NPEffect, AFDM_LRR_Var<3) #removing outliers
AFDM_NP <- rma.mv(AFDM_LRR, AFDM_LRR_Var, random = ~ 1 | Site, data=NPEffect2)
summary(AFDM_NP)

#Forestplot

library(forestplot)
tabletext <-cbind (
  c("Nitrogen (n=51)", "Phosphorus (n=45)", "Nitrogen + Phosphorus (n=41)"))

forestplot(labeltext=tabletext, 
           mean=c(0.0645, 0.2331, 0.3707), 
           lower=c(-0.0122, 0.0739, 0.2064), upper=c(0.1412, 0.3924, 0.5350),
           title="AFDM Effect Sizes by Nutrient Treatment",
           xlab="LRR",
           xticks=c(0, 0.25, 0.5, 0.75, 1),
           txt_gp=fpTxtGp(label=gpar(cex=1.25),
                          ticks=gpar(cex=1.25),
                          xlab=gpar(cex=1.25),
                          title=gpar(cex=1.5)),
           col=fpColors(box="black", lines="black", zero = "gray50"),
           cex=0.5, lineheight = "auto", boxsize=0.15, colgap=unit(6,"mm"),
           lwd.ci=1, ci.vertices=TRUE, ci.vertices.height = 0.05)

```

##Microbial Meta-Analysis AI Effect Sizes

```{r}
##Multi-level Regressions with Experiments Nested within Stream for AI (Overall Effect Sizes)

AI_N <- rma.mv(AI_LRR, AI_LRR_Var, random = ~ 1 | Site, data=NEffect)
summary(AI_N)

AI_P <- rma.mv(AI_LRR, AI_LRR_Var, random = ~ 1 | Site, data=PEffect)
summary(AI_P)

AI_NP <- rma.mv(AI_LRR, AI_LRR_Var, random = ~ 1 | Site, data=NPEffect)
summary(AI_NP)

#Forestplot

library(forestplot)
tabletext <-cbind (
  c("Nitrogen (n=51)", "Phosphorus (n=45)", "Nitrogen + Phosphorus (n=41)"))

png("AIMeta_plot.png")

forestplot(labeltext=tabletext, 
           mean=c(-0.3657, -0.1359, -0.6095), 
           lower=c(-0.5796, -0.4066, -0.9354), upper=c(-0.1518, 0.1348, -0.2836),
           title="AI Effect Sizes by Nutrient Treatment",
           xlab="LRR",
           xticks=c(-1, -0.75, -0.5, -0.25, 0, 0.25),
           txt_gp=fpTxtGp(label=gpar(cex=1.25),
                          ticks=gpar(cex=1.25),
                          xlab=gpar(cex=1.25),
                          title=gpar(cex=1.5)),
           col=fpColors(box="black", lines="black", zero = "gray50"),
           cex=0.5, lineheight = "auto", boxsize=0.15, colgap=unit(6,"mm"),
           lwd.ci=1, ci.vertices=TRUE, ci.vertices.height = 0.05)

dev.off()


tiff("AIPMeta_plot.tiff", width=12, height=12, units="cm", res=300)

forestplot(labeltext="",
           mean=c(-0.1359), 
           lower=c(-0.4066), upper=c(0.1348),
           title=expression(Effect~of~P~on~AI),
           xlab="LRR",
           xticks=c(-0.5, -0.25, 0, 0.25, 0.5),
           txt_gp=fpTxtGp(label=gpar(cex=1.5),
                          ticks=gpar(cex=1.5),
                          xlab=gpar(cex=1.5),
                          title=gpar(cex=1.75)),
           col=fpColors(box="black", lines="black", zero = "gray50"),
           cex=0.5, lineheight = "auto", boxsize=0.1, colgap=unit(6,"mm"),
           lwd.ci=1, ci.vertices=TRUE, ci.vertices.height = 0.025)

dev.off()

```

##Grazer Meta-Analysis Chlorophyll a Effect Sizes

```{r}

Grazers=read.csv("BeckHall_2018_GrazerMetaAnalysis_LogResponseRatios.csv")
NGrazers=subset(Grazers, Trt=="N")
PGrazers=subset(Grazers, Trt=="P")
NPGrazers=subset(Grazers, Trt=="NP")

##Multi-level Regressions with Experiments Nested within Study (Overall Effect Sizes)

Grazer_N <- rma.mv(LRR, LRR_Var, random = ~ 1 | Study, mods=~Grazing, data=NGrazers)
summary(Grazer_N)

Grazer_P <- rma.mv(LRR, LRR_Var, random = ~ 1 | Study, mods=~Grazing-1, data=PGrazers)
summary(Grazer_P)

Grazer_NP <- rma.mv(LRR, LRR_Var, random = ~ 1 | Study, mods=~Grazing, data=NPGrazers)
summary(Grazer_NP)



tiff("GZMeta_plot.tiff", width=12, height=12, units="cm", res=300)

tabletext <-cbind (
  c("Grazed", "Ungrazed"))

forestplot(labeltext=tabletext, 
             mean=c(-0.0166, 0.1713), 
             lower=c(-0.2406, -0.0444), upper=c(0.2074, 0.3870),
             title=expression(Effect~of~P~on~Chlorophyll~paste(italic(a))),
             xlab="LRR",
             xticks=c(-0.5, -0.25, 0, 0.25, 0.5),
             txt_gp=fpTxtGp(label=gpar(cex=1.5),
                            ticks=gpar(cex=1.5),
                            xlab=gpar(cex=1.5),
                            title=gpar(cex=1.75)),
             col=fpColors(box="black", lines="black", zero = "gray50"),
             cex=0.5, lineheight = "auto", boxsize=0.1, colgap=unit(6,"mm"),
             lwd.ci=1, ci.vertices=TRUE, ci.vertices.height = 0.025)

dev.off()

```

